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This  paper  describes  a  model  based  method  for  real  time  battery  cell  state-of-charge  (SoC)  estimation 
using  linear  parameter  varying  (LPV)  system  techniques.  For  this  class  of  methods,  the  applicable  struc¬ 
ture  is  one  in  which  the  input  to  output  dynamics  of  the  battery  can  be  described  by  a  discrete  parameter 
varying  state  variable  model  that  includes  the  SoC  as  a  state.  Within  this  context,  the  problem  of  state- 
of-charge  estimation  is  viewed  as  a  state  estimation  problem,  so  that  a  state  estimator  is  designed  using 
the  model.  An  LPV  system  technique,  combined  with  input  to  state  stability  criteria,  is  used  to  analyze 
the  stability  and  performance  of  the  estimator.  Compared  with  algorithms  available  in  the  current  litera¬ 
ture,  such  as  those  employing  an  extended  Kalman  filter  and  sliding  mode  observers,  this  method  offers 
good  performance  with  a  guarantee  of  stability,  and  possesses  user  friendly  tuning  with  low  computa¬ 
tional  complexity  for  easy  on-board  implementation.  Experimental  results  are  given  which  validate  the 
proposed  methodology. 

©  201 1  Elsevier  B.V.  All  rights  reserved. 


1.  Introduction 

As  hybrid  and  electric  vehicle  technology  advance  further,  auto¬ 
motive  manufacturers  have  adopted  lithium  ion  batteries  as  the 
electrical  energy  storage  device  of  choice  in  current  and  envisioned 
vehicles.  These  high-capacity,  high-power  batteries  provide  signif¬ 
icant  improvement  in  terms  of  energy  and  power  density  when 
compared  to  NiMH  and  lead-acid  batteries  used  in  previous  gener¬ 
ations  of  plug-in  hybrid  electric  vehicles  (PHEVs),  hybrid  electric 
vehicles  (HEVs)  and  electric  vehicles  (in  total,  these  classes  are 
referred  to  as  P/H/EVs).  However,  the  lithium  ion  chemistry  is  one 
such  that  even  mild  over-charging  or  over-discharging  can  result 
in  catastrophic  failures  or  premature  aging  [1  ].  Therefore,  the  cur¬ 
rent  strategy  to  prolong  the  life  of  the  battery  pack  is  to  use  the 
batteries  in  a  conservative  fashion.  On  the  other  hand,  utilizing  the 
batteries  to  the  full  extent  of  their  capabilities  can  provide  better 
fuel  economy,  drivability  and  reduced  cost  (for  example,  resulting 
in  a  smaller  battery  pack).  One  solution  to  simultaneously  maxi¬ 
mize  both  battery  life  and  still  garner  the  highest  utility  from  the 
batteries  is  to  have  an  effective  battery  management  system  that 
can  operate  closer  to  the  safety  limit  to  maximize  gains. 
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One  of  the  most  important  and  yet  difficult  problems  in  design¬ 
ing  an  effective  battery  management  system  is  the  estimation  of 
the  state-of-charge  (SoC)  of  the  battery  pack.  In  previous  gener¬ 
ations  of  HEVs,  the  battery  pack  is  operated  in  charge-sustaining 
mode,  usually  in  a  small  range  around  75%  (often  60-85%).  In  such 
applications,  the  margin  of  error  is  large  because  the  batteries  are 
operated  relatively  far  from  their  safety  limits.  With  the  advent  of 
PHEVs,  battery  packs  often  must  be  used  in  capacity  ranges  down 
to  around  25%  SoC  in  charge-depleting  mode  and  then  operated  in 
charge-sustaining  mode  around  25%.  Having  a  more  accurate  SoC 
estimate  means  that  even  when  the  batteries  are  near  the  limit¬ 
ing  SoC,  the  entire  hybrid  system  can  still  be  controlled  so  that 
performance  is  not  compromised  by  the  safety  systems,  which 
operate  by  curbing  battery  usage  when  voltage  limitations  are 
violated. 

Estimating  the  SoC  accurately  in  real  time  is  a  challenging  prob¬ 
lem.  Formally,  SoC  is  defined  as  the  ratio  of  available  ampere-hour 
(Ah)  to  the  total  Ah  available  when  the  battery  is  fully  charged 
(namely,  the  capacity).  When  the  current  through  the  battery  is 
measured  precisely,  the  SoC  (herein  represented  by  the  variable  z) 
can  be  calculated  via  a  Coulomb  counting  process  in  the  manner 


where  K  is  a  factor  that  that  is  inversely  proportional  to  the  capacity 
of  the  battery.  However,  because  measurement  of  electrical  current 
is  always  corrupted  by  noise,  the  integration  operation  results  in  a 
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Nomenclature 

re  disturbance  and  uncertainty  in  the  error  dynamics 
of  the  observer 

Ae  system  matrix  for  the  error  dynamics  of  the  observer 

X  slope  of  a  secant  line  on  the  OCV  to  SoC  function/ 

A T  maximum  change  of  temperature  in  one  sampling 

A{/\,  B,  C,  D}  uncertainties  in  system  matrices 
7  estimated  quantity  from  the  observer  -  i.e.  z  is  the 

estimated  SoC 

7  error  variable  between  estimated  and  true 

{A,  B,  C,  D}  state  space  system  matrices  for  the  dynamic  volt¬ 
ages  in  the  battery  model 
Ce  extended  output  system  matrix 

Cn  nominal  capacity  (Ah) 

/  open  circuit  voltage  function 

/  battery  current  (A) 

id  current  direction  (0/1  =  charge/discharge) 

Ke  overall  observer  gain 

Kx  observer  gain  for  the  dynamic  voltage  states 
Kz  observer  gain  for  SoC  state 

Lf  lower  bound  for  x 

P  positive  definite  matrix  in  the  Lyapunov  equation 

P jj  :  i  =  {0, 1},  j  =  {c,  d }  components  of  the  P  corresponding  to 
temperature  and  current  direction 
Q.  positive  definite  matrix 

T  cell  temperature 

Ts  sampling  time  (s) 

u  input  signal  to  the  model  (current) 

Uf  upper  bound  for  x 

V  Lyapunov  function  defined  for  the  observer 

Vh  hysteresis  voltage 

Voc  open  circuit  voltage 

w,  v  input  and  output  disturbances 

X  states  representing  dynamic  voltages  in  the  battery 

model 

z  state  representing  battery  SoC  (%) 


drift  of  the  estimate  away  from  the  true  value  over  time.  Addition¬ 
ally,  because  capacity  of  the  battery  degrades  over  time,  current 
integration,  which  depends  on  capacity,  will  also  become  inac¬ 
curate.  Alternatively,  one  could  measure  the  open  circuit  voltage 
(OCV)  of  the  battery  because  it  has  a  near  one-to-one  correspon¬ 
dence  with  the  SoC  and  is  independent  of  the  capacity.  However, 
battery  dynamics  dictate  that  the  OCV  can  only  be  measured  after 
the  battery  has  been  resting  for  a  significant  period  of  time.  There¬ 
fore,  the  OCV  is  not  a  useful  measurement  technique  when  the 
battery  is  under  continuous  operation.  Furthermore,  for  batteries 
whose  OCV  is  relatively  flat  with  respect  to  SoC,  this  method  can 
be  difficult  to  apply. 

Aside  from  direct  (or  indirect)  measurement,  the  SoC  can  also 
be  estimated  using  other  techniques  (for  a  brief  summary  see  [2]) 
which  improves  upon  the  shortcomings  of  the  measurement  tech¬ 
niques.  The  study  [3]  proposes  a  fuzzy  logic  method  where  the 
SoC  is  determined  using  the  frequency  response  of  the  battery. 
This,  however,  is  not  practical  in  an  on-board  setting  where  the 
controller  is  prohibited  from  injecting  a  sinusoidal  signal  into  the 
battery  pack.  Similarly,  [4]  provides  a  fuzzy  logic  current  integra¬ 
tion  method  where  the  current  integration  is  adapted  based  on 
the  operating  conditions.  Although  more  practical  than  the  pre¬ 
vious  method,  this  approach  still  suffers  from  the  drawback  of 
the  integration  process,  where  the  SoC  estimate  experiences  drift 
over  time.  Black  box  approaches  such  as  those  employing  neural 


networks  have  also  been  used  for  SoC  estimation;  for  example, 
[5]  shows  an  estimator  constructed  using  an  artificial  neural  net¬ 
work.  These  methods  can  often  produce  very  good  results  after  the 
network  is  trained  sufficiently.  However,  the  training  process  is 
laborious  (sometimes  non-convergent),  especially  if  the  estima¬ 
tor  is  expected  to  function  under  various  operating  conditions. 
Furthermore,  the  black  box  nature  makes  these  estimators  less 
intuitive. 

An  alternative  class  of  algorithms  used  for  SoC  estimation  is 
one  utilizing  model  based  techniques.  In  these  methods,  control 
theoretic  techniques  are  applied  to  a  control-oriented  model  to 
estimate  the  SoC.  Generally  speaking,  the  term  control-oriented 
usually  refers  to  a  class  of  models  that  are  of  low-order  and  have 
a  low-degree  of  nonlinearity,  yet  have  sufficient  accuracy  for  the 
subsequent  control  design  to  produce  good  results  with  the  non¬ 
linear  plant.  There  are  several  examples  of  these  algorithms  in 
the  literature.  First  is  the  extended  Kalman  filter  (EKF)  approach 
[6,7],  wherein  a  slightly  nonlinear  discrete  state  space  model  for 
the  battery  is  identified,  which  includes  the  SoC  as  a  state.  Then 
an  extended  Kalman  filter  is  applied  to  the  model  to  estimate  the 
total  state.  The  primary  drawback  of  this  approach  is  the  fact  that 
an  error  covariance  matrix  (whose  size  is  the  order  of  the  model) 
must  be  propagated  through  the  system  at  each  sampling  instance 
to  calculate  the  correction  gain.  This  results  in  a  significant  num¬ 
ber  of  calculations  that  may  not  be  suitable  for  implementation. 
In  addition,  because  the  EKF  requires  linearization  of  the  plant  at 
each  time  instance  around  an  unknown  operating  point  (namely 
the  SoC),  convergence  is  not  guaranteed. 

Another  model  based  approach  that  has  shown  success  is  the 
sliding  mode  observer  approach;  [8,9]  show  two  variations.  In 
this  approach,  a  linear  battery  model,  simplified  as  compared  to 
the  model  used  in  [6],  is  used  to  describe  the  battery  dynam¬ 
ics.  Then  a  sliding  mode  observer  is  used  to  estimate  the  SoC. 
By  allowing  the  sliding  mode  gain  to  dominate  the  plant  uncer¬ 
tainty,  this  method  is  able  to  guarantee  that  the  estimator  has 
desirable  convergence  properties.  The  estimation  results  are  good 
for  the  data  used  with  this  method.  However,  explicit  results  are 
not  given  for  regions  where  the  battery  is  operated  at  low  tem¬ 
peratures.  In  our  experience,  because  battery  parameters  (such 
as  internal  resistance)  change  significantly  when  temperature  is 
lowered,  even  when  using  a  sliding  mode  observer  whose  gains 
are  tuned  to  dominate  the  maximum  uncertainty,  the  estima¬ 
tion  errors  that  result  are  often  quite  large.  Furthermore,  because 
temperature  is  a  measured  quantity,  it  makes  sense  that  adapt¬ 
ing  the  estimator  gains  with  respect  to  temperature  can  improve 
performance. 

In  this  paper,  a  robust  SoC  estimator  is  designed  using  a  linear 
parameter  varying  (LPV)  estimator.  The  model  used  to  describe  the 
battery  dynamics  is  a  discrete  LPV  state  variable  model  which  can 
be  obtained  via  methods  discussed  in  previous  work  [10-12],  The 
parameters  on  which  the  LPV  model  operates  are  SoC  and  tem¬ 
perature;  such  parametric  variation  built  into  the  model  improves 
the  accuracy  over  various  operating  conditions  to  reduce  the  error 
in  the  estimation,  especially  at  lower  temperatures.  Because  the 
SoC  is  a  parameter  of  the  plant  and  unknown,  this  structure  treats 
the  plant  as  uncertain;  as  such,  the  resulting  estimator  must  be 
robust  to  this  uncertainty.  In  addition,  the  effect  of  the  uncer¬ 
tainty  can  be  analyzed  explicitly  using  input  to  state  stability 
criterion  [13-16].  There  are  several  advantages  to  this  design. 
First  the  computational  complexity  is  much  less  than  the  EKF 
approach  because  no  error  covariance  matrix  is  propagated  for¬ 
ward.  Second,  stability  of  the  estimator  can  be  guaranteed  while  the 
properties  of  the  estimation  error  can  be  analyzed  explicitly.  After 
the  models  are  developed  and  discussed,  the  estimator  design  is 
demonstrated  using  experimental  data  collected  from  two  different 
batteries. 
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2.  Battery  model 


In  this  section,  the  model  of  the  battery  is  described.  Because 
the  final  algorithm  is  intended  for  on-board  implementation,  the 
battery  model  is  assumed  to  be  in  a  discrete  state  variable  form 
given  as 


iwJ 


[  X[k]  \  [o  A(T,id)\  [X[k]\  + 
y[fc]  =/(T,  z)  +  CX[k]  +  D(T,  z,  id)u[fc]  + 


36C„  (u[k]  +  w[fc]) 


(2) 


where  T  is  the  temperature,  z  is  the  SoC,  id  is  the  current  direction, 
Ts  is  the  sampling  time,  C„  is  the  nominal  capacity,  X  e  n<R  is  the 
state  vector  representing  the  battery  dynamics, /is  the  open  circuit 
voltage  as  a  function  of  T and  z,  and  w  and  v  are  zero  mean  Gaussian 
noise  disturbances  that  affect  the  input  and  output,  respectively.  In 
addition,  A,  B,  C,  D  are  the  state  variable  system  matrices  of  appro¬ 
priate  dimensions  for  the  dynamic  state X.  Note  here  that  the  matrix 
A  is  assumed  to  be  independent  of  z  (the  SoC).  Previous  modeling 
results  [17]  have  shown  that  the  dependence  of  A  on  z  has  only  a 
very  minor  effect  on  the  accuracy  of  the  model.  In  addition,  because 
z  is  a  quantity  to  be  estimated,  it  must  be  considered  unknown. 
As  such,  including  this  term  does  not  improve  the  accuracy  of  the 
model. 

This  particular  battery  model  structure  represents  a  common 
approach  in  approximating  the  battery  dynamics.  It  is  based  on 
the  physical  intuition  that  the  output  voltage  of  a  battery  is  com¬ 
posed  of  an  OCV,  a  set  of  dynamic  voltages  that  are  the  results  of 
electrochemical  effects  like  charge  transferring  and  diffusion,  and  a 
voltage  resulting  from  the  internal  resistance.  The  combination  of 
the  dynamic  voltages  and  the  voltage  corresponding  to  the  internal 
resistance  is  commonly  referred  to  as  the  overvoltage  of  the  bat¬ 
tery.  An  additional  component  often  included  in  battery  models  of 
this  type  is  the  hysteresis  voltage,  which  describes  the  fact  that  the 
rested  battery  OCV  at  a  given  SoC  and  temperature  can  be  different 
depending  on  the  previous  excitation.  The  dynamic  equation  most 
commonly  used  to  describe  this  phenomenon  (see  [6,1 8]  for  exam¬ 
ple)  is  a  first  order  dynamic  equation  in  the  SoC.  For  instance,  the 
continuous  time  equation  can  be  written  as 

Vh  =  r\imz,T)-Vh),  (3) 


where  Vh  is  the  hysteresis  voltage,  r>  0  provides  the  time  constant, 
and  H  is  a  function  that  provides  the  maximum  hysteresis  voltage 
for  a  given  SoC  and  temperature.  Examining  this  equation  reveals 
that  if  the  current  input  is  zero,  then  the  equation  reduces  to 

Vh  =  0.  (4) 


The  discrete  form  of  this  equation  is 


Vh[k+\)  =  Vh[k], 


(5) 


which  is  precisely  the  same  form  as  the  dynamic  equation  for  z 
(and  consequently  for  the  OCV)  when  the  current  is  zero.  Because 
the  battery  terminal  voltage  reflects  the  combined  effect  of  the  OCV 
and  the  hysteresis,  one  cannot  distinguish  between  the  effects  of 
the  OCV  and  hysteresis  based  on  the  output  voltage  when  the  cur¬ 
rent  is  zero.  In  other  words,  if  the  hysteresis  effect  is  included  in 
the  model,  then  under  zero  current  conditions,  the  system  is  not 
observable  from  the  output.  As  such,  a  hysteresis  element  cannot 
be  included  in  a  state  estimator  design,  and  is  therefore  not  included 
in  the  remainder  of  this  paper.  Note  that  this  does  not  mean  that  the 
hysteresis  effect  must  be  ignored  altogether.  If  a  hysteresis  model 
is  available  in  real  time,  one  can  simply  remove  the  hysteresis  volt¬ 
age  computed  by  this  model  from  the  battery  voltage  measurement 
prior  to  applying  the  state  estimator.  In  so  doing,  the  influence  of 
the  hysteresis  can  be  removed  from  the  estimated  states. 


One  aspect  of  the  model  that  is  not  addressed  in  this  paper  is  the 
model  adaptation  with  respect  to  aging.  The  main  reason  for  omit¬ 
ting  this  is  the  lack  of  experimental  aging  data  at  the  present  time. 
The  problem  of  modeling  parameter  variation  as  a  result  of  aging  is 
a  topic  that  will  be  studied  as  part  of  the  future  work.  One  possible 
approach,  such  as  already  proposed  in  previous  works  involving 
the  extended  Kalman  filter  [6],  is  to  augment  the  model  states  with 
critical  parameters  such  as  internal  resistance  and  capacity.  The 
estimator  can  then  be  designed  to  estimate  these  additional  states 
as  well  as  the  SoC.  The  overall  system  is  capable  of  adapting  to  the 
aging  of  the  battery. 

The  model  described  here  can  be  identified  in  two  ways.  The 
first  method  is  to  identify  a  parameterized  equivalent  model  such 
as  done  in  [10],  In  this  case,  typically  an  equivalent  circuit  model 
consisting  of  an  OCV,  an  internal  resistance,  and  multiple  pairs  of 
parallel  RC  circuits  is  used.  The  number  of  RC  circuit  pairs  depends 
on  how  accurately  one  wants  to  model  the  battery.  Typically,  a 
second  or  third  order  model  is  sufficient  for  capturing  the  basic 
charge  transfer  and  diffusion  dynamics  in  the  frequency  range  of 
10  mHz-10  Hz.  Even  though  the  slower  diffusion  dynamics  can 
reach  as  low  as  1 0  p,Hz,  these  dynamics  are  difficult  to  model  for  this 
type  of  application,  where  the  model  structure  must  be  appropriate 
for  control  design.  For  instance,  often  utilized  elements  such  as  the 
Warburg  impedance  or  transmission  line  cannot  be  used  here  with¬ 
out  a  meaningful  time  domain  representation.  However,  the  effect 
of  these  unmodeled  dynamics  and  uncertainty  can  be  mitigated  by 
applying  control  theory.  Also  note  that  the  optimal  model  order 
may  be  slightly  different  depending  on  temperature.  However, 
having  different  model  order  for  different  temperatures  makes  esti¬ 
mator  design  very  difficult.  Therefore  it  is  assumed  that  model 
order  is  constant  for  all  temperatures.  The  model  constructed  this 
way  can  be  represented  by  a  continuous-time  differential  equation, 
which  can  be  discretized  (see  Appendix  A  for  more  information 
on  this).  Another  method  is  to  directly  identify  the  model  using  a 
black  box  identification  technique  ( one  example  can  be  seen  in  [  1 7  ] , 
where  a  subspace  method  is  used  to  perform  the  identification). 
The  model  that  results  is  very  similar  to  the  discretized  equivalent 
circuit  model.  Once  again,  lower  order  models  (such  as  second  or 
third  order)  are  typically  used  to  capture  the  basic  charge  transfer 
and  diffusion  dynamics.  For  more  information  on  the  actual  iden¬ 
tification,  interested  readers  should  refer  to  the  aforementioned 
references. 


3.  Estimator  design 


In  this  section,  the  design  of  the  SoC  estimator  is  described.  The 
basic  structure  of  this  estimator  is  a  state  observer  that  uses  the 
voltage,  temperature,  and  current  measurements  as  feedback  sig¬ 
nals.  Because  the  SoC  is  a  state  of  the  model,  if  the  states  of  the 
observer  converge  to  the  true  states,  SoC  estimation  is  achieved. 

The  generic  form  of  the  state  observer  design  is  given  as 


P  °  1  P[kll  + 

[0  A(T,id) J  [X[fc]j  ^ 


y[k]  =/(T,  2)  +  CX[k]  +  D(T,  2,  id)u[k]. 


36  Cn  1  u{k]  +  Ke(y[k]-nk]) 
_B(T,z,id)\ 


(6) 


Note  that  this  is  precisely  the  form  of  a  model  with  output  feed¬ 
back  for  correction.  The  exact  form  of  the  output  feedback  I<e  will 
be  discussed  later,  except  to  note  that  it  is  a  function  of  the  dif¬ 
ference  between  the  estimator  output  and  the  measured  voltage 
output.  Because  w  and  v  are  unknown,  they  are  not  included  in  the 
estimator.  For  estimator  coefficients  that  depend  on  the  SoC,  the 
estimated  SoC  z  is  used  in  place  of  the  true  SoC  z. 

The  output  error  y[k]  =y[k]~  y[k]  can  be  computed  as 


y[k]  =f(T,  z)  — /(T,  z)  +  CX[k]  +  A  D[k]u[k]  +  v[k], 
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where  X[k]  =X[k] -X[k]  and  AD[k]  =  D(T,z,  id)  -  D(T,z,  id).  The 
estimation  error  dynamics  can  then  be  computed  as 


rz[k+i] 

X[k+1] 


0  1 

\m' 

A(T,id)\ 

[Xlk] 

—Ts  1 

36C„ 

B(T,z,id)\ 

w[k]- 

+AD[k]u[k]  +  v[k]). 


Next,  define 

m= 


f(T[k],z[k])  — /(T[k],  z[k]) 
z[k]-z[k] 


(7) 

(8) 


For  a  given  T,  the  function/must  satisfy  a  basic  property  of  the  OCV 
function  in  that  it  is  monotonically  increasing  with  respect  to  SoC. 
Because  x  represents  the  slope  of  a  secant  line  on  the  graph  of/and 
z,  the  monotonicity  of/ implies  that  /  is  always  positive.  Because 
the  domain  of  T  and  z  are  both  compact,  3 Lf,  Uf>  0  such  that 


0  <Lf< 


xM  <  Uf, 


O) 


Vk  e  N.  Note  here  that  Lf  and  Uf  can  also  be  functions  scheduled  on 
temperature.  However,  because  the  shape  of  the  OCV  functions  is 
likely  to  be  similar  even  for  different  temperatures,  x  is  primarily  a 
function  of  the  SoC  rather  than  temperature.  Therefore  the  variation 
of  Lf  and  Uf  with  respect  to  the  temperature  will  be  small,  and  thus 
for  the  sake  of  simplicity,  Lf  and  Uf  are  taken  to  be  constants. 

Given  this, 


/(T,z)-/(T,z)  +  CX[k]  =  [X[k]C] 


X[k]  I 


In  other  words,  the  output  error  can  be  written  as  a  linear  func¬ 
tion  of  the  estimation  error,  with  the  exception  that  part  of  the 
linear  coefficient,  namely  x.  is  time  varying  and  unknown  (since 
it  involves  z  which  is  the  quantity  to  be  estimated).  This  suggests 
that  the  estimator  to  be  designed  can  be  considered  a  linear  state 
estimator  that  is  robust  to  variations  in  x-  In  other  words,  one  can 
select  Ke  as  a  linear  function  of  the  output  error  and  design  it  in  such 
a  way  that  the  estimator  converges  regardless  of  the  true  values  of 
X-  Thus,  Ke  can  be  written  as 


\Kz(T,id)  1 
'  [Kx(T,id)\ 


(10) 


where  the  components  Kz  and  Kx  could  depend  on  T  and  ij. 
The  error  dynamics  can  now  be  written  as 


(11),  even  though  Kz  and  Kx  are  assumed  to  be  functions  of  these 
scheduling  variables. 

Heuristically,  several  observations  can  be  made  immediately. 
For  linear  systems,  if  the  unforced  system  (namely  re  =  0)  has 
asymptotically  stable  dynamics,  then  with  bounded  input,  the  state 
response  will  still  be  bounded  (input  to  state  stability).  Therefore, 
the  purpose  of  Ke  is  to  make  the  unforced  system  asymptotically 
stable.  For  a  parameter  varying  system,  this  is  not  as  simple.  But 
because  the  parameter  variation  with  respect  to  T  is  very  slow, 
the  idea  of  gain  scheduling  suggests  that  placing  the  eigenvalues 
of  this  matrix  inside  the  unit  circle  will  likely  stabilize  the  sys¬ 
tem.  Along  the  diagonal  of  Ae,  there  are  two  portions:  1  -  Kzx  and 
A(T,  z,  id)  -  AA[k]  -  KXC.  If  |Xxxl  <  1  and  Kz  >  0,  then  the  eigenval¬ 
ues  of  Ae  will  be  very  close  to  1  -  KZX  and  the  eigenvalues  of  A(T, 
z,  id)  —  AA[k]  -  KXC.  Because  A  has  stable  dynamics,  the  function  of 
I(z  is  to  provide  a  greater  margin  of  stability  and  counter  effects  of 
unmodeled  dynamics.  If  this  is  accomplished,  then  the  system  will 
be  stable.  To  ensure  that  the  estimation  error  is  as  small  as  pos¬ 
sible,  re  must  be  as  small  as  possible.  With  the  exception  of  the 
last  term  of  re,  the  other  terms  are  the  result  of  model  uncertainty 
and  noise,  which  cannot  be  influenced  once  the  model  is  given.  The 
last  term  in  (13)  contains  uncertainty  and  disturbance  multiplied 
by  Ke  (in  view  of  (10)).  This  suggests  that  if  Ke  had  relatively  small 
magnitude,  greater  robustness  to  uncertainty  could  be  expected. 
However,  because  Ke  is  also  used  to  place  the  poles,  this  also  means 
that  a  faster  closed  loop  response  would  come  at  the  cost  of  higher 
sensitivity  to  uncertainty  and  disturbance. 

Theses  arguments  can  be  formalized  via  Lyapunov  theory. 
Define  a  parameter  dependent  Lyapunov  function  V  by 

V[k]=Xj[k]P(T,id)Xe[k],  (14) 

where  Pen  +  lxn+llHisa  positive  definite  matrix.  Here,  P  is 
assumed  to  depend  on  T  and  id,  even  though  if  a  choice  of  P  that 
is  independent  of  the  parameters  could  be  found,  then  it  can  and 
should  be  used. 

For  a  given  V,  we  have 

V[k+  1]  -  V[k]  =  Xer[k](4 [k]P[k+  l]Ae[k]  -P[k])Xe[k] 

+  2/^[k]P[k+l  ]Ae[k]Xe[k]+7^  [fc]P[k  +  1  ]Te[k]. 

Suppose  that  for  some  P(  ),  Q>0,  Ae[k]P[k+ l]Ae[k] -P[k]  < 
-Q,  Vk.  Note  here  that  Q  can  be  assumed  to  be  independent  of 
the  scheduling  parameters  because  if  a  parameter  dependent  Q 
exists,  then  because  the  domain  of  the  parameters  is  compact,  a 


\m+ 1]' 

[ 


where 
Ae[k]=  [ 


M1 

0 

1  \~Z[k] 1 

H° 

A(T,id) 

+Te[k], 

1-KzX 

—Kz 

C  1 

-KxX 

A(T,  id)  - 

-Xxcj 

U[k]  + 


(12) 


r'W=  [  ABlfc]  ]  “[fcl+  L|fU  W[k]-  [j]  ^Dlk]u{k\+V[k]). 

(13) 


Note  that  Ae  determines  the  dynamics  of  the  error  and  re  repre¬ 
sents  the  disturbance  or  uncertainty  that  ultimately  determines  the 
size  of  the  estimation  error.  Further  note  that  for  brevity  in  nota¬ 
tion,  the  dependence  of  Kz  and  Kx  on  T  and  id  is  not  expressed  in 


w[k]-  ([x[k]C]  +  AD[k]u[k]  +  v[k]^ 


(11) 


non-parameter  dependent  Q  can  always  be  found  to  bound  the 
parameter  dependent  Q.  It  follows  that 

V[k+  1]  -  V[k]  =  Xj [k](-Q)Xe[k]  +  2 T] [k]P[k  +  1  ]Ae[k]Xe[k] 
+7^[k]P[fc  +  l]7%[k].  (15) 

Because  the  quadratic  term  on  the  right-hand-side  of  (15) 
dominates  for  large  Xe,  3m >0  such  that  if  |Xe[fc]|  >  m,  then 
V[k+1]  — V[k]<0.  Therefore,  Xe  will  converge  to  a  neighborhood 
around  0.  The  size  of  this  neighborhood  depends  on  the  magni¬ 
tude  of  the  remaining  terms  on  the  right-hand-side  of  (15),  which 
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is  directly  dependent  on  the  magnitude  of  _Te.  This  analysis  con¬ 
firms  analytically  the  reasoning  behind  the  previous  argument  on 
the  magnitude  of  re. 

This  analysis  also  suggests  that  Ke  should  be  selected  so  that  the 
LMI  ATe[k]P[k  +  \]Ae[k]-P[k]  <  -Q  has  solutions  P{  ),  Q>0.  This 
can  be  done  using  an  LMI  solver  that  simultaneously  solves  for  all 
three  variables.  However,  such  a  process  is  generally  very  limiting 
because  Ke  and  P  must  be  in  a  specific  form  for  the  LMI  solver  to  be 
applicable.  Furthermore,  because  Ke  also  affects  the  performance 
of  the  estimator  in  terms  of  convergence  and  disturbance  rejec¬ 
tion,  the  solution  that  an  LMI  solver  finds  may  not  be  desirable.  A 
practical  solution  is  to  use  the  idea  of  gain  scheduling  to  find  Ke  and 
then  confirm  the  stability  of  the  resulting  solution  by  checking  the 
LMI  given  previously.  The  strategy  is  to  hold  T  constant  and  solve 
for  Ke  and  then  schedule  I<e  with  respect  to  T.  This  strategy  makes 
sense  because  the  temperature  process  is  slower  than  the  electrical 
process.  Therefore  the  design  problem  reduces  to  that  of  selecting 
Ke  for  various  fixed  temperatures  and  then  interpolating  the  results 
with  respect  to  temperature. 

When  designing  l<e  for  a  fixed  temperature,  the  properties  of 
X  can  be  used  to  significantly  simplify  the  design  problem.  As 
discussed  previously,  the  eigenvalues  of  Ae  should  be  very  close 
to  1  -  Kzx  and  the  eigenvalues  of  A(T,  id)  -  KXC.  The  eigenvalue 
1  -  Kzx,  which  corresponds  to  the  state  z,  determines  the  conver¬ 
gence  properties  of  the  SoC  component  of  the  observer.  Because 
X>0,  Kz  >0  is  used  to  bring  this  eigenvalue  inside  the  unit  circle. 
The  further  inside  the  unit  circle  this  eigenvalue  is,  the  better  con¬ 
vergence  property  one  has  for  the  SoC.  Because  x  is  generally  very 
small  and  Kz  should  also  be  small  to  reduce  the  effect  of  the  dis¬ 
turbance,  1  -  I<zX  is  always  positive  and  is  generally  very  close  to 
the  unit  circle.  Given  this  interpretation,  it  is  easy  to  see  that  for 
all  the  values  that  x  can  take,  X  =  Lf  results  in  the  worst  perfor¬ 
mance  for  this  design  problem.  Therefore  when  designing  Ke  for  a 
fixed  temperature,  the  problem  can  be  simplified  by  using  x=L/. 
If  the  performance  of  the  resulting  estimator  is  good  under  this 
construction,  it  will  be  even  better  for  other  values  of  X- 

The  overall  design  process  for  fixed  temperature  reduces  to  plac¬ 
ing  the  closed  loop  poles  of  the  estimator  at  a  set  of  values  that 
provide  a  good  compromise  between  convergence  speed  and  dis¬ 
turbance  rejection.  Because  disturbance  rejection  implies  smaller 
magnitude  for  the  components  of  Ke,  the  poles  selected  should  be 
very  close  to  the  original  poles.  This  is  particularly  true  for  the  state 
z  since  the  small  magnitude  of  x  requires  a  very  large  Kz  to  result  in 
a  pole  far  away  from  the  unit  circle.  Therefore  it  is  usually  advisable 
to  place  the  pole  corresponding  to  z  slightly  inside  the  unit  circle. 

After  the  pole  placement  is  done  satisfactorily,  the  LMI  condition 
resulting  from  the  Lyapunov  analysis  must  be  checked  to  verify  that 
the  closed  loop  system  is  stable.  The  LMI  in  compact  notation  can 
be  written  as 


where  AT  is  the  maximum  amount  of  temperature  change  that 
can  occur  within  one  sampling  period.  The  idea  of  (18)  is  that 
given  P(T[k  +  l],  id[k  + 1]),  P(T[k],  id[k])  must  be  inside  the  range 
P(T[k  + 1]±  AT,  {c,  d}).  However,  because  the  dependence  on  T  is 
linear,  it  is  only  necessary  to  check  the  corner  points  of  this  range. 
In  particular,  these  are  precisely  the  four  points  P(T[fc  + 1  ]  +  AT,  c), 
P(T[k  + 1  ]  +  AT,  d),  P(T[k  + 1  ]  -  AT,  c),  and  P(T[k  + 1  ]  -  AT,  d).  There¬ 
fore,  instead  of  checking  for  all  possible  P[k  + 1  ],  P[k],  it  is  sufficient 
to  only  check  these  four  possibilities.  Following  this  approach,  even 
though  the  number  of  LMIs  increases,  the  problem  becomes  much 
more  solvable. 


3.1.  A  different  look  at  SoC  estimation 

The  SoC  estimator  design  highlighted  the  two  most  important 
factors  that  influence  the  design  of  any  voltage-based  SoC  estima¬ 
tor:  the  accuracy  of  the  OCV  function  and  the  small  slope  of  the 
OCV  function.  A  natural  question  to  ask  is  whether  or  not  we  can 
estimate  the  OCV  instead  of  SoC  to  avoid  dealing  directly  with  the 
accuracy  of  the  OCV  function.  Indeed  a  similar  state  estimator  con¬ 
struction  can  be  used  to  estimate  the  OCV  based  on  the  model.  To 
do  this,  we  first  perform  a  change  of  coordinates  to  replace  the  state 
z  with  a  state  that  represents  the  OCV.  Recall  that  the  OCV  is  related 
to  the  SoC  and  temperature  via  the  function/.  Let  Voc  represent  the 
OCV.  Then 

. v 


Voc[k+  1] «  Voc[k]  -  j-(T[k],z[k])= 


(19) 


using  a  first  order  Taylor  expansion  under  the  assumption  that 
ASoC  and  AT  are  both  small  over  each  sampling  period.  Given  (19), 
the  dynamic  equation  of  the  system  can  be  written  as 


V„c[k+1] 

X[k] 


.»][ 


Voc[k] 

X[k] 


W36C„  (u[k]+w[k]) 


B(T,z,id ) 

y[k]  =  [1  C]Xe[k]  +  D(T,  z,  id)u[k]  +  v[k\. 

The  effect  of  this  nonlinear  change  of  coordinates  is  that  the  non¬ 
linear  function  /  is  moved  from  the  output  equation  to  the  input 
coefficient.  Using  this  form,  we  can  estimate  the  OCV  using  a  state 
estimator  given  as 


r  voc[k+ i]i 

[  m  \ 


i 


r  Vocik]] 

[  X[k]  J 


36C„  ulk]  +  Ke(T,id)(y{k]-y[k]) 


Ae[k]P[k+l]Ae[k]-P[k]  +  Q_  0  0 

0  -P[k]  0  <0. 

0  0  -Q 


y[k]  =  [1  C]Xe[k]  +  D{T,z,id)u[k]. 

Note  that  in  this  equation,  z  is  an  estimated  SoC  computed  using 
the  estimated  OCV  and  the  inverse  of  the  OCV  function/.  Define 


This  form  is  difficult  to  use  because  two  time  indices  occur.  There¬ 
fore  we  make  the  simplification  that  P  takes  a  linear  functional 
form: 


p(T[k],id[k])={  ;0< 

[  >0, 

Then  it  is  sufficient  1 


+  T[k]Plc  if  charging 
+  T[k]Pld  if  discharging  ' 


AeP(T,  id)Ae[k]  -  P(T,  {c,  d})  ±  ATP1{c  d}  +  Q 
0 
0 


-P(T,id)  0 
0  -Q 


*1*1 (22) 

Note  this  is  very  similar  to  (12)  with  the  exception  that  x[k]  has 
been  replaced  by  l.Just  like  Ae  in  (12),  this  matrix  determines  the 
stability  of  the  closed  loop  estimator.  Consequently,  the  same  anal¬ 
ysis  used  previously  can  be  used  to  conclude  that  if  there  exists 
symmetric  P  and  Q  such  that 

ATe[k]P[k+l]Ae[k]-P[k]  <  -Q,  (23) 

Vk,  then  the  error  dynamics  will  be  input-to-state  stable.  What 
is  different  about  this  linear  matrix  inequality  is  that  it  does  not 
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depend  on  x ■  In  fact,  if  we  assume  that  Ae  does  not  depend  on  z,  then 
there  is  no  uncertainty  in  this  matrix  inequality.  Consequently,  the 
pole  placement  approach  discussed  previously  can  be  applied  very 
easily.  After  placing  the  poles,  the  previous  matrix  inequality  can 
be  used  to  ensure  that  the  parameter  variation  of  Ae  does  not  cause 
the  system  to  go  unstable. 

The  primary  difference  between  this  estimator  and  the  previ¬ 
ous  estimator  used  for  SoC  estimation  is  that  the  feedback  is  not 
affected  by  the  quantity  Consequently,  the  pole  placement  is  not 
as  restrictive  as  is  the  case  for  the  SoC  estimation.  Thus  a  natural 
question  to  ask  is  why  not  just  estimate  the  OCV  and  then  use  the 
OCV  to  calculate  the  SoC  based  on  the  inverse  of  /?  The  answer  is 
that  while  this  is  certainly  valid,  the  problem  of  the  small  slope  of 
/  cannot  be  avoided  altogether.  When  we  calculate  the  SoC  using 
the  OCV,  the  small  errors  in  the  estimated  OCV  become  magnified 
by  the  inverse  function  of/.  In  addition,  because  the  derivative  of 
f  which  depends  on  the  estimated  SoC,  is  used  as  the  input  coeffi¬ 
cient  in  the  OCV  equation,  the  uncertainty  in  the  true  SoC  has  the 
same  effect  as  an  input  disturbance.  Therefore,  in  either  case,  the 
estimation  result  will  depend  on  the  accuracy  of  the  OCV  function 
as  well  as  the  nonlinearity  of  df/dz. 

This  alternative  estimation  scheme  does  offer  an  advantage 
when  compared  to  the  direct  SoC  estimation  scheme.  When  esti¬ 
mating  the  SoC  directly,  we  must  be  concerned  about  the  estimate 
exceeding  the  physical  limits  (over  100%  or  below  0%).  Therefore  in 
practical  implementations,  a  saturation  function  must  be  used  to 
limit  the  resulting  SoC.  Because  this  saturation  is  performed  directly 
on  a  dynamic  state,  the  stability  properties  can  be  altered.  However, 
if  the  OCV  is  estimated,  then  we  are  not  concerned  with  setting 
an  artificial  limit  on  the  resulting  estimate  because  the  measure¬ 
ment  has  the  same  units.  A  negative  voltage,  for  example,  will  only 
increase  the  feedback  error  which  would  prompt  the  estimator  to 
correct  the  estimation  accordingly. 


robustness  will  result  in  a  slower  estimator.  Hence  the  designer 
must  decide  which  attribute  is  more  important.  Furthermore,  this 
conflict  also  couples  with  the  problems  that  arise  when  Lu  is  small. 
In  such  cases,  a  faster  convergence  requirement  can  significantly 
degrade  the  tracking  robustness.  Design  examples  are  provided  in 
the  next  section  to  illustrate  these  issues  more  clearly. 

4.  Design  examples 

In  this  section,  the  design  algorithms  described  previously  are 
illustrated  using  data  collected  from  two  separate  batteries.  Due 
to  chemistry  differences,  the  estimator  behavior  is  very  different 
in  these  two  cases.  The  differences  highlight  the  design  tradeoff 
described  in  the  previous  section. 


4.1.  EIC  battery 


In  this  section,  a  constant  temperature  design  example  is 
provided  using  data  collected  from  an  EIG  battery  that  has  a 
capacity  of  20  Ah  and  nominal  voltage  of  3.6  V.  This  battery  uses 
a  Li[NiCoMn]C>2  based  cathode  and  a  graphite-based  anode.  The 
model  of  the  battery  is  identified  using  the  data  collected  from  the 
battery  operating  under  an  asymmetrical  step  profile.  This  model 
includes  both  an  OCV  map  as  a  function  of  the  SoC  and  the  dynamic 
system  matrices  required  by  (2 ).  The  OCV  is  measured  by  repeatedly 
discharging  the  battery  by  10%  and  then  resting  for  1  h  to  record 
the  OCV  for  that  SoC.  The  SoC  of  the  battery  for  the  current  pro¬ 
file  is  computed  by  integrating  a  post-processed  current  (noise  and 
sensor  offsets  are  removed). 

The  A  and  C  matrices  for  the  model  equation  (2)  are  given  as 


(24) 


3.2.  Comments 

In  the  above,  the  design  of  the  SoC  estimator  (as  well  as  an  OCV 
estimator)  is  described  in  mathematical  detail.  Some  important 
points  that  might  otherwise  be  lost  among  the  equations  should  be 
re-emphasized.  First,  the  analysis  shows  that  for  any  voltage-based 
SoC  estimator,  the  quantity  that  ultimately  determines  how  well 
the  estimator  will  work  is  Lu,  the  smallest  slope  of  the  OCV  func¬ 
tion.  Larger  Lu  will  make  the  SoC  more  sensitive  to  voltage  variation, 
which  allows  us  to  obtain  more  from  the  voltage  measurement 
signal.  Second,  there  is  always  a  conflict  between  convergence 
speed  and  tracking  performance.  Increased  convergence  speed  will 
inevitably  decrease  tracking  robustness,  while  increased  tracking 


The  SoC  equation  for  the  system  (accounting  for  the  capacity  and 
the  sampling  time)  is  given  by 

z[k  +  l]=z[k]-2^gU[fc].  (25) 

Figs.  2-4  show  the  input  coefficients  bi,  f>2  and  the  internal  resis¬ 
tance  D  as  functions  of  SoC  and  current  direction,  respectively.  Fig.  5 
plots  the  open  circuit  voltage  as  a  function  of  the  SoC. 

Given  these  coefficients,  a  SoC  estimator  can  be  designed  using 
the  method  described  in  the  previous  section.  First  evaluate  the 
derivative  of  the  OCV  function  /  shown  in  Fig.  6.  As  we  can  see 
from  this  plot,  this  derivative  is  bounded  between  0.004  V/%  and 
0.01 1  V/%.  Therefore  Ke  can  be  designed  for  /  =  0.004.  In  other 
words,  Ce  =  [0.004,  C], 


Rl 


Rn 


Fig.  1.  Equiv 
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The  next  step  is  to  select  the  pole  locations.  Aside  from  the  inte¬ 
grator  pole  in  the  state  z  (which  is  on  the  unit  circle),  the  other  two 
poles  are  quite  fast  already.  For  example,  the  pole  at  0.9877  will 
cause  an  initial  condition  on  x,  to  reduce  to  less  than  3%  of  its  value 
in  300  samples  (150s),  which  is  quite  fast  for  vehicle  applications. 
Therefore  there  is  no  need  to  select  poles  for  these  two  states  to 
be  much  faster  than  their  current  values.  Moreover,  the  penalty  of 
having  faster  poles  is  higher  gains,  which  causes  the  system  to  be 
more  reactive  to  measurement  errors  and  noise.  Consequently,  we 
select  [0.96,  0.72]  as  the  pole  locations. 

The  main  question  for  the  design  is  where  to  put  the  pole  cor¬ 
responding  to  the  SoC.  Clearly,  the  pole  must  be  moved  inside  the 
unit  circle.  The  effect  of  the  pole  location  can  be  seen  using  a  sim¬ 
ulation.  First,  artificial  Gaussian  zero  mean  random  noise  is  added 
to  the  current  and  voltage  measurement.  Then  additional  pulses 
of  various  magnitude  and  time  duration  are  added  to  the  current 
measurement.  This  simulates  non-zero-mean  measurement  errors 
that  can  mislead  a  simple  current  integrator.  Fig.  7  shows  an  exam¬ 
ple  of  the  input  disturbance.  Note  that  the  size  of  the  disturbance 
is  quite  large  so  that  its  effect  is  very  significant  on  the  system. 


Derivative  of  the  OCV  function 


SoC  [%] 

Fig.  6.  Derivative  of  the  OCV  function. 
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Fig.  7.  Disturbance  added  to  the  input  current. 


Placing  a  pole  at  0.99  results  in  kz  =  8.8740.  The  result  of  this 
feedback  gain  is  that  the  estimator  converges  very  quickly  from  an 
initial  error,  evident  in  Fig.  8  where  the  estimator  is  used  to  correct 
the  SoC  trajectory  that  was  initialized  incorrectly.  After  only  20 
samples  ( 1 0  s),  the  estimator  is  within  5%  of  the  true  SoC.  However, 
the  cost  of  this  is  that  the  input  and  measurement  noise  can  often 
cause  the  estimation  error  to  be  quite  large,  which  is  seen  in  Fig.  9. 

Next  a  much  slower  pole  of  0.9991  is  used,  which  results  in 
kz  =  0.8945.  As  expected,  the  convergence  of  the  estimator  is  much 
slower.  As  shown  in  Fig.  10,  the  estimator  converges  from  a  very 
large  initial  error  to  within  5%  of  the  true  SoC  in  1500  samples,  or 
approximately  750  s.  However,  as  seen  in  Fig.  11,  the  disturbance 
now  has  a  much  smaller  effect  on  the  estimation  error. 

Note  that  stable  experimental  performance  does  not  necessarily 
indicate  stability  under  all  operating  conditions.  Therefore,  it  must 
be  confirmed  that  the  designs  are  stable  for  all  other  values  of  /.  To 
do  this,  P  and  Q must  be  found  such  that  (18)  holds  for  all  x  e  [A/,  Uf]. 


estimation  error 


estimation  error 


Fig.  8. 
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To  this  end,  the  Matlab  LMI  toolbox  is  used  to  find  feasible  solutions 
for  P  and  Q,  For  the  pole  at  0.99,  the  following  P  and  Q  result: 

r  0.1  -7.5  2.1  1  r  0.0041  -0.6407  -0.1400  1 

P=  -7.5  5512.8  -69.3  <2=  -0.6407  125.8615  23.2645  . 

L  2.1  -69.3  2080.1  J  |_  -0.1400  23.2645  480.5619  J 

Similarly,  for  the  pole  at  0.9991,  we  have 

[  7.5  136.2  9.2  1  I"  0.0163  -0.9209  -0.1554 1 

136.2  7020.4  -91.1  <2  =  -0.9209  140.2092  21.0865  . 

9.2  -91.1  2555.4J  [_  —0.1554  21.0865  604.1492J 

The  fact  that  P  and  Q,  exist  (in  both  cases)  concludes  that  the  esti¬ 
mator  is  stable  using  either  set  of  gains. 

Given  that  both  designs  are  stable,  it  is  up  to  the  user  to  select 
the  design  that  best  satisfies  the  design  objectives.  For  vehicle  appli¬ 
cations,  we  note  that  very  large  SoC  error  is  rare.  Therefore,  there 
is  very  little  need  for  fast  convergence,  and  robustness  to  sensor 
noise  or  offset  is  much  more  important.  Consequently,  the  robust 
performance  offered  by  the  second  design,  albeit  with  slower  con¬ 
vergence,  is  a  better  solution. 

4.2.  A) 23  26650  cells 

In  this  section,  a  multi-temperature  design  example  is  pro¬ 
vided  using  data  gathered  from  an  A123  26650  lithium  ion 
iron-phosphate  battery,  which  has  a  LiFeP04-based  cathode  and 
graphite-based  anode.  This  particular  battery  has  a  capacity  of 
2.1 5  Ah  and  a  nominal  voltage  of  3.3  V.  For  this  battery,  a  model 
is  identified  for  operating  temperatures  between  — 5°C  and  45  °C 
and  for  SoC  between  10%  and  90%.  Once  again  the  model  contains 
a  measured  OCV-SoC  function.  One  simplification  that  is  made  is 
that  the  OCV  is  selected  to  be  independent  of  temperature.  This 
is  due  to  lack  of  data  points  in  lower  SoC  regions  for  lower  tem¬ 
perature  data.  This  simplification  necessarily  introduces  estimation 
errors  for  lower  temperature  data.  Nevertheless,  when  a  more  pre¬ 
cise  OCV  function  becomes  available,  it  can  simply  be  used  in  place 
of  this  simplified  OCV  function  to  yield  better  results. 

The  datasets  used  with  this  battery  are  asymmetrical  step  pro¬ 
files,  containing  current  steps  of  2  C,  4  C,  and  6  C  that  are  designed 
to  allow  the  battery  SoC  to  traverse  between  10  and  90%.  The  same 
profile  is  used  for  all  temperatures  in  the  range  described  pre¬ 
viously.  Note  that  the  range  limits  on  temperature  and  SoC  are 
dictated  by  the  availability  of  experimental  data.  For  example,  on 
a  battery  cycler,  it  is  difficult  to  design  a  profile  that  reaches  0%  or 
100%  effectively.  Often  because  the  current  control  is  not  exact,  a 
profile  that  tries  to  reach  0%  or  100%  can  cause  the  undervoltage 
or  overvoltage  safety  systems  to  stop  the  experiment  before  com¬ 
pletion.  Given  that  the  battery  is  never  operated  outside  of  the  10 
to  90%  range  in  vehicle  operation,  we  limit  the  SoC  to  be  within 
this  range.  Nevertheless,  if  the  dataset  contains  operation  in  other 
ranges,  the  methodology  discussed  here  applies  equally. 

Once  again,  the  first  task  is  to  evaluate  the  derivative  of  the  open 
circuit  voltage  function  with  respect  to  SoC;  the  plot  of  the  deriva¬ 
tive  is  given  in  Fig.  12.  First  note  that  the  lowest  point  on  the  plot 
occurs  in  two  sections:  between  35%  and  55%  SoC  and  between 
70%  and  85%  SoC,  where  the  lowest  value  reached  is  0.0004  V/%. 
Compared  with  the  same  plot  for  the  EIG  battery  seen  in  Fig.  6, 
this  number  is  almost  exactly  one  order  of  magnitude  smaller. 
Consequently,  this  suggests  that  SoC  estimation  for  this  battery  is 
inherently  a  much  more  difficult  problem  than  with  the  EIG  bat¬ 
tery.  Furthermore,  a  large  difference  in  magnitude  exists  between 
the  section  where  the  lowest  points  occur  and  the  remaining  sec¬ 
tions.  This  suggests  that  the  SoC  estimation  will  be  more  accurate 
when  the  battery  is  operating  outside  of  these  two  zones.  This  is 
somewhat  welcome  news  because  if  the  battery  starts  to  approach 
the  low  and  high  SoC  regions  where  the  d//dz  is  large,  the  estimator 


slope  of  the  open  circuit  voltage  function  f 


should  be  able  to  estimate  the  SoC  more  accurately,  thus  allow¬ 
ing  the  energy  management  system  information  to  avoid  taking 
the  battery  further  into  the  low  and  high  SoC  regions.  Flowever, 
this  also  has  an  undesirable  aspect  because  the  regions  between 
35%  and  55%  SoC  and  between  70%  and  85%  are  precisely  where 
vehicles  operate  in  charge-sustaining  modes.  The  region  between 
70%  and  85%  is  the  region  in  which  FIEVs  typically  operate  while 
the  region  around  35%  is  frequented  by  PFIEVs  after  they  enter  the 
charge-sustaining  mode.  Consequently,  this  is  not  the  ideal  battery 
for  applying  the  SoC  estimator  design.  Some  of  the  results  later  will 
reflect  this  analysis. 

For  the  purpose  of  the  estimator  design,  Lf  from  (9)  can  be 
selected  as  0.0004.  The  two  eigenvalues  of  the  A  matrix  as  func¬ 
tions  of  the  temperature  are  given  in  Table  1 .  This  table  is  obtained 
by  evaluating  the  eigenvalues  of  the  A  matrices  in  the  identified 
model  at  each  temperature  in  the  table.  We  select  the  poles  of  the 
closed  loop  matrix  as  [0.99,  0.91  ],  which  are  smaller  than  the  min¬ 
imum  values  of  the  two  poles  as  functions  of  the  temperature.  The 
reason  for  this  choice  is  that  we  require  the  closed  loop  system  to 
respond  faster  than  the  open  loop  system,  but  only  slightly  because 
the  state  z  will  be  relatively  slow.  The  location  for  the  pole  corre¬ 
sponding  to  z  is  selected  as  0.996.  While  a  faster  pole  would  provide 
faster  convergence,  because  the  If  is  so  small,  a  very  large  feedback 
gain  would  be  required  to  achieve  faster  performance.  Therefore 
we  opted  for  a  slower  pole  and  thus  a  smaller  feedback  gain. 

To  see  the  performance  of  the  estimator  with  the  choice  of  the 
poles  noted  above,  we  utilize  the  step  profile  data  that  was  used 
for  modeling.  For  the  25  °C  step  profile  data,  Fig.  13  shows  the 
estimated  SoC  compared  with  the  SoC  calculated  from  current  inte¬ 
gration  (noted  as  “measured”).  Fig.  14  shows  the  estimation  error 
and  the  input  disturbance.  From  these  two  figures,  we  can  see  that 
because  this  battery  has  a  much  flatter  OCV  function  than  the  EIG 
battery  example  we  have  seen  previously,  the  estimator  is  much 

Table  1 

Eigenvalues  of  the  A  matrix  as  functions  of  T. 

Temperature  (°C)  Pole  1  Pole  2 

45  0.9920  0.9657 

35  0.9916  0.9477 

25  0.9906  0.9399 

15  0.9916  0.9367 

5  0.9921  0.9255 

0  0.9972  0.9229 

-5  0.9971  0.9126 
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measured  and  estimated  SoC  (25°C  data)  measured  and  estimated  SoC  (45°C  data) 


Fig.  13.  SoC  estimated  using  the  25  °C  data  compared  with  the  measured  data. 

Fig.  15.  SoC  estimated  using  the  45  °C  data  compared  with  the  measured  data. 


more  susceptible  to  disturbance  and  modeling  errors.  However,  the 
overall  performance  at  room  temperature  is  still  very  good  as  the 
estimator  provides  accuracy  with  an  error  of  less  than  5%,  except 
under  offset  disturbance. 

In  Figs.  1 5  and  1 6,  the  same  estimator  is  used  with  the  45  °  C  step 
profile  data.  Because  the  modeling  error  is  even  smaller  at  45  °C,  the 
estimator  performance  at  this  temperature  is  even  better  than  with 
the  25  °C  data. 

Figs.  17  and  18  show  the  estimator  performance  over  the  5°C 
step  profile  data.  In  this  example,  we  see  that  the  performance  is 
worse  than  the  room  temperature  case.  There  are  two  reasons  for 
this.  First,  the  modeling  error  at  lower  temperatures  is  much  higher 
than  that  at  higher  temperatures.  Consequently,  the  magnitude  of 
the  uncertainties  that  influence  the  size  of  the  dynamic  error  is 
much  higher.  Second,  the  OCV  function  for  lower  temperatures  is 
different  from  the  OCV  function  at  higher  temperatures.  Because 
df/dz  is  very  small  for  this  battery,  what  might  be  considered  small 
differences  in  OCV  can  influence  the  SoC  estimation  significantly. 
To  investigate  this  further,  we  further  moved  the  closed  loop  pole  of 
z  out  toward  the  unit  circle  to  decrease  the  effect  of  modeling  error 
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Fig.  14.  SoC  estimation  error  and  input  disturbance  with  25  °C  data. 


Fig.  16.  SoC  estimation  error  and  input  disturbance  with  45  °C  data. 
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Fig.  18.  SoC  estimation  error  and  input  disturbance  with  5  °C  data. 


Fig.  20.  SoC  estimation  error  and  input  disturbance  with  a  slower  pole. 


on  the  estimator  performance.  As  Figs.  19  and  20  show,  the  slower 
pole  dramatically  reduced  the  convergence  speed.  However,  the 
tracking  performance  after  convergence  improved  dramatically. 

An  example  will  serve  to  illustrate  the  operation  and  perfor¬ 
mance  of  the  OCV  estimator.  The  25  °C  data  is  used  for  this  purpose. 
As  we  can  see  from  Fig.  23,  the  estimated  OCV  converges  to  the 
OCV  calculated  based  on  the  measured  SoC,  and  tracks  well  after 
a  period  of  convergence.  Therefore  the  OCV  estimator  discussed  in 
the  previous  section  is  a  valid  alternative  to  the  SoC  estimator. 

The  results  here  suggest  that  for  this  particular  battery,  the  OCV 
as  a  function  of  temperature  should  also  be  described.  Otherwise, 
the  accuracy  of  the  lower  temperature  estimation  is  significantly 
affected.  Due  to  absence  of  a  better  OCV  function,  for  further  illus¬ 
tration  a  multiple  temperature  example  is  provided  that  uses  model 
based  simulation  data.  In  this  example,  a  current  profile  that  fea¬ 
tures  simultaneous  temperature  and  SoC  variation  is  used  to  excite 
the  model;  subsequently,  the  estimator  is  used  to  estimate  the 
SoC  using  the  model  output  as  inputs.  As  always,  random  signals 
are  added  to  both  the  input  and  output  to  simulate  model  uncer¬ 
tainty  and  noise.  As  Figs.  21  and  22  show,  despite  the  temperature 


samples  [k] 

Fig.  21.  SoC  using  the  multi-temperature  design. 
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measured  and  estimated  OCV  (25°C  data) 


variation  and  the  uncertainty,  approximately  the  same  perfor¬ 
mance  in  terms  of  estimation  error  results.  After  the  estimator 
converges,  the  overall  estimation  error  is  kept  well  within  5%  SoC. 


5.  Conclusion 

In  this  paper,  a  SoC  estimator  is  designed  using  linear  param¬ 
eter  varying  system  techniques.  The  advantages  of  this  design  are 
threefold.  First,  this  estimator  does  not  require  iterative  calculation 
of  the  feedback  gain  such  as  in  designs  using  the  extended  Kalman 
filter.  Because  of  this,  the  computational  burden  is  relatively  small 
so  that  implementation  is  more  feasible.  Secondly,  the  stability  of 
the  estimator  can  be  confirmed  analytically  based  on  linear  matrix 
inequalities  formulated  using  input-state  stability  criterion.  Lastly, 
the  performance  of  the  estimator,  in  terms  of  convergence  and 
tracking,  can  be  tuned  very  easily  depending  on  user  requirements. 

A  fundamental  fact  emerging  form  this  work  is  that  the  most 
important  component  of  a  voltage  based  SoC  estimator  designed 
using  state  estimation  is  the  open  circuit  voltage  of  the  battery. 
This  is  the  only  component  aside  from  direct  current  integration 
that  provides  the  user  with  a  reliable  measure  of  the  SoC.  For  bat¬ 
teries  such  as  the  EIG  battery  whose  OCV  has  relatively  steep  slope, 
a  small  error  can  be  tolerated  since  the  voltage  changes  significantly 
enough  when  the  battery  SoC  changes.  But  for  batteries  such  as  the 
A123  (whose  characteristics  are  common  to  all  lithium-ion  iron- 
phosphate  chemistry),  the  inherent  flatness  of  the  OCV  makes  it 
critical  to  have  as  good  of  an  OCV  model  as  possible.  Consequently, 
an  important  continuation  of  this  work  is  to  improve  the  OCV 
model,  especially  as  a  function  of  temperature.  Such  a  model  can 
be  obtained  with  a  well  designed  measurement  experiment.  Even 
though  such  an  experiment  may  be  tedious  or  time  consuming,  its 
impact  on  the  accuracy  of  the  estimator  cannot  be  underestimated. 
Furthermore,  the  data  from  this  experiment  would  lead  to  formu¬ 
lation  of  a  hysteresis  model,  which  combined  with  an  accurate  OCV 
model  would  provide  even  better  estimates. 


Appendix  A.  Equivalent  circuit  model  discretization 

The  equivalent  circuit  model  used  in  [10]  is  provided  in  Fig.  1. 
This  model  comprised  of  an  open  circuit  voltage,  an  internal 
resistance,  and  multiple  R/C  circuits  that  represent  the  battery 


(A.1) 


dynamics.  The  dynamics  of  each  parallel  RC  circuit  ai 
by  a  first  order  differential  equation  of  the  form 

^  =  -AM  +AiB,I, 
at 

where  A,  and  B,  represent  the  input  coefficient  and  the  time  con¬ 
stant,  respectively.  Note  that  A,  and  Bj  can  be  functions  of  the 
parameters  as  well,  but  for  the  purpose  of  discretization,  it  is 
assumed  that  the  parameters  do  not  change  within  each  sampling 
instance. 

Because  the  OCV  is  a  static  function  of  the  SoC,  discretizing  the 
OCV  dynamics  can  be  done  by  discretizing  the  dynamics  of  the  SoC. 
The  continuous-time  dynamics  of  SoC  denoted  by  z  are  given  by 


i  3600 C„1, 

where  C„  is  the  nominal  capacity  of  the  battery. 

The  overall  battery  terminal  voltage  is  given  by  the  si 
the  components 


(A.2) 


vhatt  =  Voc  -  RqI  - 


Given  a  selected  sampling  period  Ts  and  a  continuous  linear 
system  of  the  form 

x  =  ax  +  bu,  (A.4) 

where  x,ueVt  and  the  coefficients  a,  b  and  the  input  u  are  assumed 
to  be  constant  over  each  sampling  period,  the  discrete  equivalent 
of  the  system  is  given  by 

x[k+  1]  =  e~aTsx[k]  +  |(1  -e~aT*)u[k].  (A.5) 

Using  this  fact,  Eqs.  (A.l )  and  (3)  can  be  discretized  as 

Vi[k+1]  =  aiV[k]  +  biI[k],  (A.6) 

Vh[k  +  1]  =  Y{I[k])Vh[k]  +  K(I[kmz,  T),  (A.7) 

where 

a,-  =  e-W, 

bj  =  B,(l  -e-AA), 

y(/[k])  =  e-r|IIfcH, 

K(/[k])  =  1  -  e-rM. 


Eq.  (A.2)  simply  represents  an  integrator,  and  can  therefore  be 
discretized  as 


z[k+l]=z[k]~  - 


(A.8) 


The  complete  model  can  then  be  written  in  state  variable  form  as 


'  z[k+  1]  ■ 

V,  [k  +  1  ] 

"1  0  . 

0  a j  . 

..  0  0 

..  0  0 

■  z[k\  ' 
Vi  [fc] 

V„[k+1] 

0  0  . 

. .  a„  0 

Vnjk] 

_Vh[k+l]_ 

0  0  . 

..  0  y(/[k])_ 

vhW. 

Ts 

3600Cn 

b\ 


(A.9) 


where  the  output  equation  (battery  terminal  voltage)  is  given  by 


Vb[k]  =  Voc(z[k])  -  R0[k]I[k]  -  £Vj[k]  -  V„[k],  (A.10) 
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After  discretization,  this  system  becomes  essentially  indistin¬ 
guishable  from  the  discrete  model  form  (2). 
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